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Abstract 

We present a new algorithm for reconstructing an unknown source in Thermoacoustic and 
Photoacoustic Tomography based on the recent advances in understanding the theoretical nature 
of the problem. We work with variable sound speeds that might be also discontinuous across 
some surface. The latter problem arises in brain imaging. The new algorithm is based on an 
explicit formula in the form of a Neumann series. We present numerical examples with non- 
trapping, trapping and piecewise smooth speeds, as well as examples with data on a part of the 
boundary. These numerical examples demonstrate the robust performance of the new algorithm. 

1 Introduction 

Thermoacoustic (TAT) and Photoacoustic (PAT) Tomography are emerging medical imaging modal- 
ities [32, 30]. These are hybrid medical imaging methods that combine the high resolution of 
acoustic waves with the large contrast of optical waves. TAT and PAT have been developed to 
overcome the limitations of both conventional ultrasound and microwave imaging. The physical 
principle underlying TAT and PAT is the photoacoustic effect which can be roughly described as 
follows. A short impulse of electromagnetic microwaves or light is sent through a patient's body. 
The tissue heats up slightly and the heat expansion generates weak acoustic waves. These waves 
are measured away from the patient's body and one tries to recover the acoustic source that gives 
us information about the rate of absorption at each point in the body, thus creating an image. 
One of the potential applications for TAT is early breast cancer detection. The American Cancer 
Society reports that breast cancer is the second overall leading cause of death among women in the 
United States. The mortality rate from breast cancer has declined in recent years due to progress 
in both early detection and more effective treatment. Better detection techniques, however, are 
still needed. The significance of TAT and PAT is that they yield images of high electromagnetic 
contrast at high ultrasonic resolution in relatively large volumes of biological tissues [32, 30]. 
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A first step in TAT and PAT is to reconstruct the amount of deposited energy from time- 
dependent boundary measurement of acoustic signals. To start with, we describe here the widely 
accepted mathematical model of TAT and PAT [32, 30, 6]. Let O C M n be an open set with 
a smooth strictly convex boundary. Let c(x) > be the sound speed that is either smooth or 
piecewise smooth. Assume that c(x) = 1 outside O. The acoustic pressure u(t, x) then solves the 
wave equation 

' (d| - c 2 A)u = in (0,T) x E n , 

u\t=o = f, (1) 
d t u\ t =o = 0, 

where T > is fixed, and f(x) is a source that we want to recover supported in fl. The measure- 
ments are modeled by the operator 

A / : = u \[o,T]xdn- (2) 

Given A/, the problem is to reconstruct the unknown / that is related to the absorbing properties 
of the body at any point. 

There have been significant progresses in both mathematical theories and medical applications; 
see [1, 6, 7, 9, 10, 12, 13, 16, 17, 19, 23, 27, 28, 32, 35, 36, 14] and references therein. Theoretically, 
one is interested in uniqueness and stability of the solution for the inverse problem; numerically, 
one is interested in designing efficient numerical algorithms to recover the solution of the inverse 
problem. Naturally, the above two aspects have been well studied in the case of the sound speed 
being constant. In fact, if the sound speed is constant and the observation surface dQ, is of some 
special geometry, such as planar, spherical or cylindrical surface, there are explicit closed-form 
inversion formulas; see [6, 31, 9, 10, 5] and references therein. In practice there are many cases 
when the constant sound speed model is inaccurate [32, 14, 34, 18]. For instance in breast imaging, 
the different components of the breast, such as the glandular tissues, stromal tissues, cancerous 
tissues and other fatty issues, have different acoustic properties. The variations between their 
acoustic speeds can be as great as 10% [14]. 

To tackle variable sound speeds in TAT and PAT, the time-reversal method has been suggested 
in [6] and used in [33, 8, 12, 13]. However, only when the dimension is odd and the sound speed 
is constant, the time-reversal method gives an exact reconstruction for T large enough by the 
Huygens principle. When the dimension is even or the sound speed is not constant, the time- 
reversal method only yields an approximate reconstruction as T > 1. A natural question arises 
immediately: is there any exact reconstruction formula which is able to handle both a variable 
sound speed and an irregular observation geometry, for a fixed T? It turns out that the work 
[27] summarized below does provide such a formula and this formula is the foundation for our 
algorithmic development. Some results about variable sound speeds are contained in the references 
above but a more complete analysis in this case has been given in [27]. This analysis includes 
if-and-only-if conditions for uniqueness and stability, including the case with observations on a part 
of the boundary. An explicit recovery formula of a type of convergent Neumann series is derived in 
[27] when A/ is known on the whole <9fi and T is greater than the stability threshold. This formula 
is the base for our numerical reconstruction, and we also apply it if the speed is trapping. 

There is another situation related to a variable discontinuous sound speed which arises in brain 
imaging [35, 20]. The skull has a discontinuous sound speed which is piecewise smooth with jump- 
discontinuities across the boundary of the skull. A typical situation is that the speed in the skull is 
about twice as large as that in the soft tissue; see e.g., [35, 36]. Such speeds change drastically the 
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way that singularities propagate; see Figure 1. In [28], the second and the third author studied this 
model and proved that the Neumann series expansion works as well, provided that all singularities 
issued from supp / have a path of non-diffractive segments that reaches d£l until time T. Based on 
this work, we develop an efficient numerical algorithm as well to handle the situation of variable 
discontinuous sound speeds in PAT and TAT. 

The rest of the paper is organized as follows. In Section 2, we provide some preliminary mate- 
rials for understanding theories summarized in the following sections. In Section 3 we summarize 
the theoretical results in [27] and explain consequences of the theory. In Section 4 we summarize 
the theoretical results in [27] for the case of partial data. In Section 5 we summarize the theoretical 
results in [28] when the sound speed is discontinuous. In Section 6 we make comparison between 
smooth and nonsmooth sound speeds in terms of uniqueness and stability. In Section 7 we summa- 
rize the iterative algorithm for constructing the Neumann series so that the inversion formula in 
[27] can be implemented. In Section 8 we give an extensive examples to demonstrate the robustness 
of the new algorithm for PAT and TAT. 

2 Preliminaries 

Assume for now that c > is smooth. The speed c defines a Riemannian metric c~ 2 dx 2 . For any 
piecewise smooth curve [a, b] 4 7 G P, the length of c in that metric is given by 



The so defined length is independent of the parameterization of 7. The distance function dist(x, y) 
is then defined as the infimum of the lengths of all such curves connecting x and y. 

For any (x, 9) £ W 1 x 5 n_1 we denote by 7^,0 (i) the unit speed (i.e., I7I = £(7)) geodesies issued 
at x in the direction 9. 

Recall that the energy of u(t, x) in a domain U is given by 



where u(t) = u(t, •). The energy of any Cauchy data (/, g) for equation (1) is given by the same 
integral with V x u = V x f and u t = g. In particular, the energy of (/, 0) in U is given by the square 
of the Dirichlet norm 



We always assume below that / € H d(^), where the latter is the Hilbert space defined by the norm 
above. We will denote by || • |[ the norm in Hd{Q), and in the same way we denote the operator 
norm in that space. 

There are two main geometric quantities that are crucial for the results below. First we set 



Let Xi < 00 be the supremum of the lengths of all maximal geodesies lying in Cl. Clearly, T < Xi 
but while the first number is always finite, the second one can be infinite. It can be shown actually 
that 






Tq := max{dist(x, dQ); x € £1}. 



(3) 



T < Xi/2. 



(4) 
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One of the main ingredients of our approach in [27, 28] is understanding the microlocal nature 
of the problem. We recall the definition of a wave front set of a function, or more generally, a 
distribution; see [11]. The definition is based on the known property of the Fourier transform: one 
can tell whether a compactly supported function / is smooth by looking at the decay of the Fourier 
transform /(0 as |£| -> oo: / G Cg° if and only if < C N (l+\£\)- N for any N. The idea behind 

the wave front set is to localize this near a fixed xq and in a conic neighborhood of a fixed £o 7^ 0. 
Conic neighborhoods are defined as open conic sets, i.e., sets of the type T = {r9; r > 0, 9 G V}, 
where V is an open subset of S" 1-1 . We say that (xo,£o) WF(/), £o 0, if there exist <p G Co° 
with (j)(xo) 7^ and a conical neighborhood T of £o> so that 

\M{0\<c N {i + \£\)- N , veer,viv. 

If (x,£) G WF(/), we say that (x,£) is a singularity of /, or that / is singular at (x,£). Since 
singularities are defined by conic sets, we can restrict £ to unit vectors. For example, the Dirac 
Delta function <5(x) has wave front at x = and all directions, i.e., WF(<5) = {(0, £); £ ^ 0}, and a 
piecewise smooth function / that has a jump across some smooth surface S (and nowhere else) is 
singular at all points of S in (co)normal directions, i.e., WF(/) = {(x,£); x G S 1 , 7^ £ ± S" at x}. 

The propagation of singularities for the wave equation (1) can be described as follows. If 
(x,9) G WF(/), then at time t, both (7x,e(*),7x,e(*)) and (7x,e(-*),7x,e(-*)) are in WF(u(t,-)), 
where u is the solution of (1). This is due to the fact that the symbol of the wave operator has two 
sound speeds, ±c(x)|£|, and that the initial velocity on (1) is zero, therefore each singularity splits 
into two equal parts starting to propagate in opposite directions. While this is a classical result in 
the linear PDE theory, we refer to [27] for more details in this specific case. 

3 Smooth speed and data on the whole <9Q 

We will describe below the theoretical results in [27]. 

3.1 Uniqueness 

If T S> 1, A/ recovers / uniquely. We have the following sharp result based on the unique contin- 
uation theorem by Tataru [29]. 

Theorem 1. Let Af = 0. Then f(x) = for dist(x, dCl) < T. Moreover, f can be arbitrary in the 
set dist(x, <9J1) > T , if the latter set is non-empty. 

Corollary 1. A is injective on Hrj(^i) if and only ifT> To. 

We refer to [27] for proofs. 

3.2 Stability 

We showed in [27] that Af recovers / in a stable way if each singularity i.e., each element of 

the wave front set WF(/), reaches dQ for time t (positive or negative) such that \t\ < T. In other 
words, if functions / are a priori supported in a fixed compact Kcfi, then we have the following 
equivalent statement: 

For any (x, 9) G K. x S n ~ l , the unit speed geodesic through (x, 9) at t = 

(5) 

reaches oQ at time \t\ < T. 
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Moreover, the following condition is sufficient regardless of the choice of KL: 

Ti/2 < T. (6) 

This condition is equivalent to (5) if /C = Q. 

Furthermore, to formulate a condition equivalent to (5) by taking K, into account, we define 
T\ = Ti(/C) as above which is explicitly related to geodesies that pass through JC. 

If Ti = oo, then the sound speed c is called trapping (in SI). In this case, there is no stability 
regardless of the choice of T. 

We summarize the above discussion into the following. 

Theorem 2. Let K, C SI be compact. 

(a) Let T > T\{K,)/2. Then there exists a constant C > so that 

ll/H < C||A/||„ 1([0)T]xan) . 

(b) Let T < Ti(/C)/2. Then for any C > 0, s± and S2, there is f £ C°° supported in K, so that 

> C\\Af\\ H s 2 ^ 0T ] xdn y 

In other words, T > Ti(/C)/2 is a sufficient and necessary condition for stability for functions / 
supported in K, up to replacing the < sign by <. 

Visible and Invisible Singularities. The condition (5) can be explained in the following 
way. As a general principle, for a stable recovery in such a linear inverse problem, we need to detect 
all singularities. We refer to [26] for more details. As explained above, each singularity starts to 
travel in a positive and a negative direction because ut = at t = 0, so it can leave two traces on 
dSl. It is enough to detect one of them for stability. On the other hand, if we can detect both, we 
can expect better numerical results. Condition (5) then says that all singularities are visible at the 
boundary. 

We want to emphasize that T\ can be much larger that diam(fi) := max{dist(x, y); (x,y) € 
dSl x dSl}. If there are no conjugate points in SI, then those two quantities coincide. If there are 
conjugate points however, this is not necessarily so. If there are closed geodesies, then T\ = co, 
while diam(fi) is finite. In this case, although we always have uniqueness for some T > 1, we never 
have stability. Finally, we notice that while To and diam(fi) can always be estimated analytically 
or numerically, T\ is much harder to estimate or even to tell whether it is finite or not. In any case, 

diam(ft) < Ti. (7) 

If Ixfi does not hit the boundary at time \t\ < T, we call (x, 6) an invisible (possible) singularity. 
It is easy to show that if one such pair exists, then there is a non-empty open set of invisible 
singularities. 

3.3 Reconstruction 

The reconstruction method in [27] is based on the following ideas. If we knew the Cauchy data 
(u, ut) on {T} x SI, we could just solve a mixed problem like the one below with that Cauchy data 
on t = T and boundary data given by Af. Although we do not know (u, u t ) on {T} x Q, we do 
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know the boundary values of u on dQ for t = T, i.e., u on {T} x dQ. Assuming that [/, 0] is in the 
energy space, i.e., / £ Hr)(^l), we can only say that Uf(T,-) is in L 2 (Q), and its boundary values 
might not be well defined. Now from all possible functions with prescribed boundary values on 
{T} x 0, we choose the one that minimizes the energy norm || • \\H D (n)- By the Dirichlet principle, 
it is given by the harmonic extension of u(T, -)\dn- 

Consequently, given h defined on [0,T] x d£l which eventually will be replaced by A/, we first 
solve the elliptic boundary value problem 

A<p = 0, <P\ dn = h(T r ), (8) 

and we introduce the notation Pq for the Poisson operator of harmonic extension: PqIi(T, •) := </>. 
In fact, the equation is c 2 Au = but c 2 cancels out. Then we perform the modified back projection: 

( (3 t 2 -c 2 A)v = in (0,T) x n, 

v \[o,T]xdn = h, , . 

v\ t =T = Pnh(T,-), w 

V t \t=T = 0. 

Note that the initial data at t = T satisfy compatibility conditions of first order (no jump at 
{T} x dfl). Then we define the following left pseudo-inverse of A 

Ah:=v(0,-) inn. (10) 

The operator A is not an actual inverse (unless n is odd, c is constant, and T is greater than the 
diameter) and we have 

AA = I - K, 

where K is an "error" operator. We showed in [27] that 

11^/11^(0) < \\f\\ HD{ ciy V/ e H D (n), (11) 

for any smooth speed, trapping or not, and for any time T > 0. If T > To, the inequality is strict, 
see [27]. To show that K is a contraction requires some assumptions, however. 

Theorem 3. 

(a) Let c be non-trapping, and let T > T\/2. Then AA = I — K, where \\K\\ HD (Q)-^H D (n) < 1- 
In particular, I— K is invertible on Hjj(Q,), and the inverse thermoacoustic problem has an explicit 
solution of the form 

oo 

f=Y,K m Ah, h:=Af. (12) 

(b) Let T > T\. Then in addition to the conclusions above, K is compact in Hrj(^l). 
For a proof, see [27, Theorem 1] and [28, Remark 2.2]. 

In other words, under condition (6), we have not only stability but an explicit solution in a 
form of a convergent Neumann series as well. On the other hand, if (6) fails, i.e., if T < Ti/2, there 
is no stability by Theorem 2. This does not mean that the series (12) would not converge in this 
case. If it does, it will converge to / for T > To. Indeed, then it is easy to see that the limit g 
solves (I — K)(g — f) = 0, and since (11) is strict in this case, / = g. 
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Based on this theorem and its proof, we can expect good convergence when T > T±, and the 
first term would already be a good approximation of the high frequency part of /. 

If Ti/2 < T < T±, we can expect the first term in the series to recover only a fraction of the 
high frequency part of /, and the successful terms to improve this gradually; the series (12) would 
still converge but that convergence would be slower. 

If T < T±, then \\K\\ = 1. Indeed, if we assume that \\K\\ < 1, we would get uniform convergence 
of the Neumann series, and stability; on the other hand, there is no stability in this case. 

3.4 Summary: Dependence on T 

(i) T < T 

Af does not recover / uniquely; sec Theorem 1. Then ||i^|| = 1, and for any / supported in 
the inaccessible region, Kf = f. 

(ii) T <T < Ti/2 

This can happen only if there is a strict inequality in (4). Then we have uniqueness but not 
stability. In this case, \\K\\ = 1, \\Kf\\ < \\f\\, and we do not know if the Neumann series 
(12) converges. If it does, it converges to /. 

(iii) Ti/2 < T < Ti 

This assumes that O is non-trapping for c. The Neumann series (12) converges exponentially 
but maybe not as fast as in the next case. There is stability, and || < 1. 

(iv) Ti < T 

This also assumes that Q is non-trapping for c. The Neumann series (12) converges exponen- 
tially. There is stability, \\K\\ < 1, and K is compact. 

3.5 Recovery of singularities 

In cases (iii) and (iv) above, we recover explicitly the whole /, including its singularities. If the 
goal is to recover only the singularities of / (the wave front set WF(/)) with less computation, 
then one can perform the classical back-projection (time reversal) as follows; see [33, 8, 12, 13, 19]. 
Assuming the non-trapping condition, T\ < oo, let % £ Co°(^) be such that x(i) = 1 for t £ [0, Ti]. 
Set 

Rf = A x Af. (13) 

In general, Rf is not close to / unless T — >■ oo; see [12]. 

In case (iv), R is a parametrix of infinite order; see [27]. Therefore, it recovers correctly all 
singularities, including jumps across smooth surfaces — it will recover correctly the location and 
the size of the jump. 

In case (iii), R is elliptic but not a parametrix itself; see [27, Theorem 3]. Then Rf will have 
the singularities at the right places but the amplitudes will be in general between 1/2 and 1. 

In cases (i) and (ii), only singularities "close enough to the boundary" will be recovered with 
amplitudes between 1/2 and 1. 

Finally, if the speed is trapping in Q, i.e., if Ti = oo, then Rf recovers the visible singularities 
up to time T by choosing \ appropriately. 

These comments are just another way to formulate [27, Theorem 3]. 
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3.6 Comparison with the Time Reversal Method 

Let Ti < oo first, i.e., assume that c is non-trapping. Let the time reversal approximation be 
defined as in (13) with x = 1 in a neighborhood of [0, Ti]. Then R is a parametrix of infinite 
order, i.e., Rf = f — Qf, where Qf € C°° for any /. On the other hand, ||Q|| is not necessarily 
small for any fixed T. Assuming that \ is properly chosen, as T — > oo, that norm gets smaller at 
a rate dictated by the local energy decay for the wave equation: ||Q|| = 0(t 1 ~ n ) for n even and 
IIQH = 0(e~ Ct ) for n odd; see [12]. Therefore, for T > 1 so that ||Q|| < 1, one can write 

RAf = (I - Q)f, 

and one solve this equation by Neumann series as above. However, it is not straightforward to 
impose sharp conditions on \ an d T to guarantee ||Q|| < 1. On the other hand, for the method 
that we propose, that condition is T > Ti and it is sharp for a stable inversion. Also, the pro- 
posed method minimizes the norm of the "error" operator, and when both Neumann expansions 
converge uniformly, the one in Theorem 3 will converge faster in the uniform topology. Numerical 
experiments not shown here confirm that. 

If Ti = oo (c is trapping), then the error in the time reversal method decays like 0(1/ log T) if 
/ € Hq(Q), and the error decays even slower if / € He>(£1) only. The first term of the Neumann 
series inversion has an error no less than that, and numerically the error improves with a few more 
terms. We do not know whether it converges or not, however. 

4 Smooth speed and data on a part of dft 

Let r C dQ be a relatively open set of dQ, and assume that we only have data available on [0, T] x T. 
We will suppose that / is supported in some compact Kcfl. 

4.1 Uniqueness 

As in (3), set 

T := T (7C, r) = max{dist(x, T); x G £}. (14) 
Theorem 1 has the following analog in this case. 

Theorem 4. Let Af = on [0,T] x F. If T > T , then f = 0. If T < T , then f = on 
K. fl {x; dist(x, T) < T} and can be arbitrary in the complement of this set in K.. 

The proof of this theorem is not easy. It combines Tataru's uniqueness theorem with arguments 
that first appeared in [6] in the case of constant speed and were extended later in [27] to variable 
speeds. 

As a corollary, A/| [0 ,t] 

x r determines uniquely / if and only if T > To. 

4.2 Stability 

The following condition guarantees that we can detect all singularities originating from K, as sin- 
gularities of our data: 

For any (x,9) € K, x S" 1-1 , the unit speed geodesic through (x,9) at t = 
reaches T at time \t\ < T. 
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Let Ti/2 < oo be the maximum of all such times. The following condition then is equivalent to 
(15) 

Ti/2 < T. (16) 

Theorem 5. Let K, C £1 be compact. 

(a) Let T > Ti/2. Then there exists a constant C > so that 

([o,T]xr)- 

(b) Let T < Ti/2. Then for any C > 0, si and S2, there is f G C°° supported in K so that 

WfWnn > C||A/||^ 2 ([o,T]xr)- 

In other words, T > Ti/2 is a sufficient and necessary condition for stability for functions / 
supported in K. up to replacing the < sign by <. 

4.3 Reconstruction 

An explicit formula of the type (12) is not available in this case but one can show that the recovery 
is reduced to a Fredholm equation if (16) holds. 

Let x ^ C°°([0,T] x r) be a cutoff function supported in [0, T] x T so that x = 1 on a slightly 
smaller set [0, T'] x T' that still satisfies the condition. Assume also that < % < 1. Then we know 
xA/. Apply the time reversal operator A to that to get AxAf. Since x(">T) = 0, the harmonic 
extension is zero, so this is the classical time reversal. Let K be the "error operator" 

K = I - A X A. 

To find /, we need to solve 

(I - K)f = h, h:= A X Ah (17) 

with h determined by the data. By [27, Theorem 3], A%A is a pseudo-differential operator that is 
elliptic under the stability condition (15). Its principal symbol is given by 

1 1 

where t±(x,£) is the positive/negative time of the (unit speed) geodesic j X £ through (x, £) to reach 
Oft. Then K is also a pseudo-differential operator with principal symbol 

<? P {K) = 1 - i X (7,, ? (r+(x,0)) - ^X(7^(r-(x,0))- (18) 

Since < % < 1, and by the stability condition, at least one of the x terms above is equal to 1/2. 
Therefore, 

< a p (K) < I 

and in the set where x is n °t nor 1 (that is relatively small in applications, but not too small 
to keep \dtx\ + l^xXl under control), a p (K) is either or 1/2. By the Garding inequality, one can 
express K in the form 

K = K\ + K2, \\Ki\\ < K2 compact, 
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and if T ^ 0, one can show that actually ||-Ki|| = 1/2. This does not allow us to claim that 
the Neumann series (12) converges. However, by [22], since the distance from K to the compact 
operators is less than one, for a fixed /, 

oo 

J2 K m A x Af 

m=0 

converges if and only if K m AxAf — > as m — s> oo. It is trivial to see that the limit g solves 

A x Ag = A x Af. (20) 

While we know that %A is injective for T > To, this is not true for A, so we cannot draw the 
conclusion that g = f ■ The original equation x^-f = h with h in the range of A is Fredholm with 
a trivial kernel, but the equation (20) for g is Fredholm which may have a non-trivial kernel. 
On the other hand, the partial sums 

N 

g N := K m A X Af (21) 

m=0 

used in the reconstruction below, recover singularities of / asymptotically as N — > oo. This follows 
from the standard pseudo-differential calculus. 

For the purpose of the numerical reconstruction, we take x to be a function of x G dO, only, 
independent of t. In particular, we no longer have x = near t = T. Then we define A as before 
but this time the harmonic extension of xA/(T, •) is not trivial in general. Then we use partial 
sums as in (21). 

If the stability condition (16) is not satisfied, then inverting %A is not equivalent to solving an 
Fredholm equation anymore and it is unstable. One can show that K is the sum of an operator 
with norm not exceeding 1 and a compact one. Of course, this does not guarantee convergence of 
the Neumann series and the partial sums recover asymptotically the visible singularities only. In 
the numerical example below, we check the norm of each successive term in the partial sum (21) 
and stop when that norm starts to increase. 

5 Discontinuous sound speed; modeling brain imaging 

We will briefly review the results in [28]. Let c(x) be piece-wise smooth with a non-zero jump across 
one or several smooth closed non- intersecting surfaces that we call S in Jl. We still assume that 
c = 1 outside The formulation of the problem is still the same and classical solutions u(t, x) are 
assumed to be in C 1 , and this implies the following transmission conditions across S where c(x) 
jumps: the limits of u(t, x) and its normal derivative match as x approaches S from either side. 

5.1 Propagation of singularities 

It is well known that this case is quite different from the previous one from the viewpoint of 
propagation of singularities. Let us see what happens when a "ray" (a geodesic) approaches S from 
one of its sides, that we call an interior one. Let c_ be the limit of c on S from inside, and let c+ 
be that from outside. 



(19) 
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If c_ > c+, this ray splits into two parts when hitting S. One of them reflects according to 
the usual laws of reflection and goes back to the interior. Another one goes into the exterior and 
refracts by changing its angle with respect to the boundary. The incoming and the outgoing angles 
u± with the normal satisfy SnelFs Law: 

$^ = C -. (22) 
sin a + c + 

Assume now that c_ < c+. Then there is a critical angle < ao < vr with the normal at any point 
so that if a_ < ao, there are still a reflected and a transmitted (refracted) ray as above satisfying 
Snell's law. If a_ > ao, then there is no refracted ray, while the reflected one still exists. This is 
known as a full internal reflection. This critical angle can be read off from the Snell Law: it is the 
value of a_ that forces sina+ to be equal to 1, and therefore to be greater than 1 when a_ > ao- 
Therefore, ao = sin _1 (c_/c+). Propagation of singularities when a_ = ao is more delicate and will 
not be analyzed here. 

After splitting or not into two rays, each segment may split into two, etc. Assuming that we 
never get rays tangent to S, reflection and refraction occur in the same way as above. Figure 1 
below shows a possible time evolution of a single singularity, both for positive and negative time. 
The sound speed in the "skull" is higher than that on either side. 

We define the distance function dist(x,y) as the infimum of the length (in the metric c~ 2 dx) of 
all piecewise smooth curves connecting x and y that intersect S transversely. 



5.2 Uniqueness 

Theorem 1 and Corollary 1 still hold in this case; see [28, Proposition 5.1] and its proof. The 
definition of Tq in this case is the same. 



5.3 Stability. Visible and invisible singularities 

The general principle that we should be able to detect each singularity for stable recovery still 
applies. However, in the case of full internal reflection, there might be singularities that never 
reach dfl. For example, let O be the ball B(0,R) (centered at with radius R), and let S be the 
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sphere \x\ = R\, R\ < R. Assume also that inside the sphere S, the speed is a constant less than 
outside S, where it is also a constant. Then all rays coming from the interior of S, hitting S at 
an angle smaller than the critical one, will reflect and give rise to no transmitted rays. By the 
rotational symmetry, those reflected rays will hit S again at the same angle, fully reflect, etc. Such 
singularities will be invisible, no matter what T we choose. This is similar to the trapping situation 
in the case of a smooth sound speed. Of course, we may also have rays that are trapped but never 
reach the interface S. Therefore the singularities that are certain to be visible up to time T consist 
of the following set: 

U ={(x,6) e(n\S)x S™- 1 ; there is a path of the "geodesic" issued from either 
(x, 9) or (x, —9) at t = never tangent to S, that is outside A at time t = T}. 

If we restrict our attention to / a priori supported in some compact K C O \ 5, then the stability 
condition can be formulated as follows: 

V(x,0) G K x S' 11 ' 1 ; there is a path of the "geodesic" issued from either 

(x, 9) or (x, —9) at t = never tangent to S, that is outside £l at time t = T. 

The reconstruction formula below shows that this condition is sufficient not only for stability but 
also for an explicit convergent formula of the type (12). 

Similar to the smooth case, we can define T±/2 (instead of defining T\ directly) to be infimum 
of all T satisfying (24) when K = Cl. Then (6) remains a necessary and sufficient condition (up to 
replacing < by <) for stability in this case as well. 

5.4 Reconstruction 

In [28], we analyzed also the energy at high frequencies that is carried by the reflected and the 
transmitted rays. We showed that if none of the incident, reflected and transmitted rays are tangent 
to S, then a positive fraction of the energy reflects and transmits at high frequencies. Therefore, 
even though under the condition (24) the corresponding singularity is visible at d£l, only a fraction 
of its would be measured. This fraction could be quite small, depending on the speed, the number 
of reflections and refractions before the ray hits dQ, and the angles there. So we should not expect 
the first term in the series (12), even if it converges, to be a good approximation even at high 
frequencies. The theorem below, proved in [28], shows that an analog of (12) still converges but it 
needs to be modified first. 

The needed modification is connected to the following problem. When applying A to Af with 
supp / C /C, we get AAf that may be supported everywhere in fl. The successful terms in (12) 
would then force AA to be applied to the result. To "restrict" AAf to /C, we cannot just restrict in 
the usual sense, because this may take the results out of the energy space (no zero trace on dtC). 
For this reason, we project orthogonally on Hd(1C). It turns out that the orthogonal projection is 
given by Yljcf := f — Pjc(f\dK.), where Pic is the Poisson operator of harmonic extension from <9/C 
to /C defined in (8). 

Theorem 6. Let K, satisfy (24). Then II^Ai = I - K in H D {K), with \\K\\ Hd{ jq < 1. In 
particular, I — K is invertible on Ho(IC), and A restricted to Hd{K,) has an explicit left inverse of 
the form 

oo 

f=J2 K m U K Ah, h = Af. (25) 

m=0 
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Note that the theorem does not say how to reconstruct / when supp / is not in a set /C satisfying 
(24). Estimate (11) still holds but we do not know if I — A" can be inverted by a Neumann 
series. In the numerical examples, however, we work with / supported everywhere in Q, and then 

= n n = i. 

6 Comparison between the smooth and the non-smooth case 

If the sound speed c is smooth, then there is always uniqueness for large enough T because To < oo; 
however, there is stability and the Neumann series converges to a solution only when c is non- 
trapping and T > T\. In the trapping case, one still has stability and an explicit solution of the 
type in Theorem 3 under the a priori assumptions that supp / C JC and all singularities with a base 
point x in K, are visible. The assumption supp/ C K, means that we need to know / outside /C; 
then we can subtract that part from / and apply the reconstruction procedure. 

When c is of jump type, we may still have rays that are trapped even if they never reach S. On 
the other hand, S may split some rays. Under the stability condition (24) which is equivalent to (6) 
with the modified definition of Ti, we still have stability and an explicit solution for supp / C K.. 
We still need to know / outside /C for a full reconstruction. The essential difference is that even 
though all singularities except those (of measure zero) contributing to tangent rays on S will be 
detected in the best possible case, T\ < oo, they might be detected with a loss of energy. In the 
example presented in Figure 1, the singularity that exits at the top has been split twice before 
that, and each such event takes away a fraction of the energy. In contrast, in the case of a smooth 
sound speed, such a singularity will exit for t > with half of the energy of / while the other 
half will exist for negative time; as such, for example, we will detect both in case (iv). Going back 
to the non-smooth case, the Neumann series (25) will gradually restore the right strength of the 
singularity, and the whole /, actually, but we can expect this to be more sensitive to noise and 
computational errors. 

7 Algorithmic formulations 

The main issue in implementing the above formulation is how to compute the operator K. In the 
case of data available on the whole d£l, by definition, K = Id — AK. Thus, given a function tp, 

Ktp = ip - A Aip. 

Thus we need an efficient computational algorithm to carry out the actions of time-reversal operator 
A and the measurement operator A. Since the measurement operator A can be simulated by forward 
modeling, we will detail the implementation of forward modeling first. 

7.1 Complex Scaling/Perfectly matched layers (PML) for the acoustic wave 
equation 

Since the forward model equation (1) is formulated as a pure initial value problem, we have to 
truncate the computational domain to be finite. The truncated domain has to be large enough 
to enclose the domain where the measurements are taken along its boundary. On the other hand, 
we have to impose some artificial boundary conditions on the boundary of the truncated domain. 
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To perform accurate long-time simulation, we will adopt the complex scaling method (see [25] and 
the references therein), also known as perfectly matched layer (PML) as an absorbing boundary 
condition for acoustic waves [2, 4, 21]. This is done so that waves will not be reflected into the 
computational domain. 

To derive the PML for the equation, we follow [21]. We rewrite the 2nd-order equation as a first- 
order system by introducing the two-component velocity vector v(x,y,t) = (v x (x,y,t),v y (x,y,t)): 

% = - v "- < 26 > 

at = " cV ' v - (27) 

By using the complex coordinate stretching [2, 4, 21], we have the following equations: 

dv v du . . 

■m+w = -aij> (28) 

= -c 2 ^, (29) 

where r\ = x,y, u = + u^ y \ and u v represents a loss in the PML and is zero in a regular 
non-PML region. 

In our simulation, we will take the following uj v : [0, 1] — >■ R + : 



0, s <G [a, 1 — a]; 

^) 2 , se[l-a,l], 



where a > and b > are some appropriate constants. Depending on the computational domain, 
this function will be rescaled accordingly. 

Equations (28) and (29) are discretized by a staggered finite-difference scheme, and the details 
are omitted; see [21]. 

7.2 To and the time reversal operator A 

To compute the critical time To, we first solve the following eikonal equation, 

c\Vf\ = 1, (30) 
T|r = 0, (31) 

by using the fast sweeping finite-difference scheme as designed in [37, 15, 24], where T can be the 
whole boundary d£l or a part of 9f2. Then To = max{T}. 

The Poisson equation defining the harmonic extension is solved by a V-cycle multigrid method 
[3], which converges with the optimal rate independent of the mesh size. 

The time-reversal wave equation (9) with terminal and boundary conditions is solved by a 
standard second-order finite-difference time-domain scheme. 
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8 Numerical results: smooth speed, data on the whole dQ 

In this section we show two dimensional numerical examples to validate our algorithms. 

Our computational domain is taken to be [— 1.5, 1.5] 2 , and Q = [—1.28, 1.28] . We work with 
three sound speeds: 

ci(x,y) = 1.0 + 0.2sin(27rx) + 0.1 cos(2vry), (32) 

c ^ x ^ = i + 9( a tf ^2) + ex P (- 90 (* 2 + y 2 )) ~ °' 4ex P (~10(3V^ + V 2 ) ~ 2) 2 ) , (33) 
c 3 (x,y) = 1.25 + sin(2.07rx)cos(2vry). (34) 

See Figure 2, with some geodesies plotted. We also use a smooth cutoff for a smooth transition of 
those speeds to 1 near <9f2 and outside Q. 



The velocity The velocily The velocily 




Figure 2: Sound speed models, (a): the variable non-trapping sound speed c\. (b): the variable 
radial trapping sound speed C2- (c): the variable trapping speed C3. 



In all examples below, we use the abbreviation NS for the Neumann series, and TR for time 
reversal. We take a finite number, k + 1, of terms in the Neumann series expansion by stopping 
when the error gets below 5% in the stable cases, and somewhat higher in the non-stable ones. 
Taking more terms, especially in the cases where we have stability, improves the error several times 
but we do not show those examples. Next, we only compute and comment on the 1? error since all 
test images model non H 1 functions. We also show a diagram of the distance function, dist(x, dQ). 



8.1 Example 1: The Shepp-Logan phantom 
8.1.1 Non-trapping speed c\, Figures 3 5 

The sound speed is given by equation (32); see Figure 2 (a). The original image as the source 
function / is shown in Figure 3(b). Numerically, we estimate To to be Tq ~ 1.1767. A rough 
estimate of T% is 3 < T\ < 4. We added four bright disks (/ = 1 there) to the classical Shepp-Logan 
phantom. 

Figure 3: T = 2T w 2.3535. The time T is slightly above the stability threshold Ti/2 but 
below T\. The error of the NS solution is 6.63% with k = 8 (9 terms) of the series vs. error 37.76% 
for the TR one. Since T < T±, the TR solution does not recover the correct size of the jumps - 
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they are recovered with amplitudes ranging from 1/2 to 1, and for many of them, it is just 1/2; this 
is clear from the slice diagrams. In contrast, the NS solution has the right amplitudes and would 
improve with more terms. 

Figure 4: T = 4Tq. The time T is doubled, and it is greater than T\. The error of the 
NS solution is 4.99% with k = 8 (9 terms) of the series vs. error 7.07% for the TR one. The 
TR reconstruction gets better as expected. While the singularities are recovered correctly, the 
low frequency part in the TR reconstruction is still not well recovered as evidenced in the slice 
diagrams. For example, in Figure 4(e) (g) for the TR reconstruction there are oscillations in some 
flat regions; in contrast, in Figure 4(f) (h) for the NS reconstruction those oscillations are gone in 
the same flat regions. 

Figure 5: T = 4To with 10% noise. The time T is doubled, and it is greater than T\. The 
error of the NS solution is 7.07% with k = 8 (9 terms) of the series vs. error 9.71% for the TR one. 
The TR reconstruction gets better as expected. While the singularities are recovered correctly, the 
low frequency part in the TR reconstruction is still not well recovered as evidenced in the slice 
diagrams. For example, in Figure 5(e) (g) for the TR reconstruction there are oscillations in some 
flat regions; in contrast, in Figure 5(f) (h) for the NS reconstruction those oscillations are gone in 
the same flat regions. 

8.1.2 Trapping sound speed C3, Figure 6 

The sound speed is C3 given by equation (34) and Figure 2 (c). Although it is hard to show that 
the speed is actually trapping, one can easily construct numerically geodesies of length at least 20; 
therefore, based on this numerical evidence, T\ is at least 20 and possibly 00. The times T that 
we choose are much smaller than 20/2 so there is a considerable set of singularities that are not 
detected. 

We have To ~ 1.25. In Figure 6, we present reconstructions with T = 4T> The phantom in the 
middle of the NS image (error 18.11%, k = 8) can be separated from the background, while in the 
TR one (error 39.56%) — much harder. The invisible singularities in this case are distributed in a 
chaotic way and what appears as noise in both images is due to the fact that there are singularities 
there that cannot be resolved — even if the actual image is not singular there! 

8.2 Example 2: Zebras, Figures 7-11 

8.2.1 Non-trapping speed c\ 

The sound speed is given by equation (32) and it is visualized in Figure 2(c). Here / now is 
represented by the zebras image that has more complex structure. 

Figure 7: T = 4Tq. Both methods give a good reconstruction, and with k = 2 only, the NS 
one has error 4.6% that is about 1/2 the of the TR one. The error can be improved significantly 
with more terms. 

Figure 8: T = 4T with 10% noise. Here, k = 8, with an error 6.25%, about 1/2 the of the 
TR one. 

8.2.2 Trapping sound speed C3 

The sound speed is given by equation (34). As we indicated above, To ~ 1.23, T\ > 30, and the 
geodesic flow is somewhat chaotic with invisible singularities that have also chaotic distribution. 



Qian, Stefanov, Uhlmann and Zhao 



Thermo- and Photo-acoustic Tomography 



17 




Figure 3: Example 1 with the non-trapping speed c\. Case 2: T = 4To. (a): the boundary distance 
map. (b): the exact initial condition, (c): the time reversal solution, (d): the Neumann series 
solution, (e): x-slices of the time reversal solution (continuous line) and the exact solution (a 
dashed line), (f): x-slices of the Neumann series solution (continuous line) and the exact solution 
(a dashed line), (g): y-slices of the time reversal solution (continuous line) and the exact solution (a 
dashed line) . (h) : y-slices of the Neumann series solution (continuous line) and the exact solution 
(a dashed line). 
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Boundary distance map 



The exact initial condition 



(a) 



(c) 




(e) 



(f) 





(g) 



(h) 
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Figure 4: Example 1 with the non-trapping speed ci. Case 2: T = 4To. (a): the boundary distance 
map. (b): the exact initial condition, (c): the time reversal solution, (d): the Neumann series 
solution, (e): x-slices of the time reversal solution (continuous line) and the exact solution (a 
dashed line), (f): x-slices of the Neumann series solution (continuous line) and the exact solution 
(a dashed line), (g): y-slices of the time reversal solution (continuous line) and the exact solution (a 
dashed line) . (h) : y-slices of the Neumann series solution (continuous line) and the exact solution 
(a dashed line). 
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Boundary distance map 



The exact initial condition 




( )• 



(a) 



-1 -0.5 



0.5 1 



(b) 



-1 -0.5 



The time reversal solution 



(c) 




The Neumann series solution 



(d) 




-1 -0.5 




(e) 



(f) 











-__JL . ... 




1 

0.8 
0.6 
0.4 




















0.2 












„1 


-0.5 0.5 1 


/l \ -1 -0.5 0.5 1 

(h) 



(g) 



Figure 5: Example 1 with the non-trapping speed c\. Case 2: T = ATq. (a): the boundary distance 
map. (b): the exact initial condition, (c): the time reversal solution, (d): the Neumann series 
solution, (e): x-slices of the time reversal solution (continuous line) and the exact solution (a 
dashed line), (f): x-slices of the Neumann series solution (continuous line) and the exact solution 
(a dashed line), (g): y-slices of the time reversal solution (continuous line) and the exact solution (a 
dashed line) . (h) : y-slices of the Neumann series solution (continuous line) and the exact solution 
(a dashed line). 
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Boundary distance map 



The exact initial condition 
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The time reversal solution 
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Figure 6: Example 1 with the trapping speed C3. T = ATq. (a): the boundary distance map. (b): 
the exact initial condition, (c): the time reversal solution, (d): the Neumann series solution, (e): 
x-slices of the time reversal solution ("-") and the exact solution ("o"). (f): x-slices of the Neumann 
series solution ("-") and the exact solution ("o"). (g): y-slices of the time reversal solution ("-") 
and the exact solution ("o"). (h): y-slices of the Neumann series solution ("-") and the exact 
solution ("o"). 
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Boundary distance map The exact initial condition 




/ \ -1 -0.5 0.5 1 /l \ -1 -0-5 0.5 1 

(g) * (h) 

Figure 7: Example 2 with the non-trapping speed c\. Case 1: T = 4Tq. (a): the boundary distance 
map. (b): the exact initial condition, (c): the time reversal solution, (d): the Neumann series 
solution, (e): x-slices of the time reversal solution ("-") and the exact solution ("o"). (f): x-slices 
of the Neumann series solution ("-") and the exact solution ("o"). (g): y-slices of the time reversal 
solution ("-") and the exact solution ("o"). (h): y-slices of the Neumann series solution ("-") and 
the exact solution ( "o" ) . 
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Boundary distance map The exact initial condition 
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Figure 8: Example 2 with the non-trapping speed c\. Case 2: T = 4To with 10% noise, (a): the 
boundary distance map. (b): the exact initial condition, (c): the time reversal solution, (d): the 
Neumann series solution, (e): x-slices of the time reversal solution ("-") and the exact solution 
("o"). (f): x-slices of the Neumann series solution ("-") and the exact solution ("o"). (g): y-slices 
of the time reversal solution ( "-" ) and the exact solution ( "o" ) . (h) : y-slices of the Neumann series 
solution ( "-" ) and the exact solution ( "o" ) . 
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Boundary distance map The exact initial condition 




Figure 9: Example 2 with the trapping speed C3. Case 1: T = 4To. (a): the boundary distance 
map. (b): the exact initial condition, (c): the time reversal solution, (d): the Neumann series 
solution. 

Figure 9: T = 4Tq. The NS reconstruction is much cleaner, error 9.64% with k = 20 that is a 
bit less than 1/2 of the TR error. 

Figure 10: T = 4To with noise. The noise increases slightly the error in both cases. 

8.2.3 Trapping speed C2 

Figure 11: The sound speed is given by equation (33) and we estimate To to be To ~ 2.1547. It 
is trapping with the circles of radii approximately 0.23 and 0.67 being stable trapped rays. The 
invisible singularities are much more structured now and are in some neighborhoods of those two 
circles. As a result, jumps across radial and close to radial lines near those circles are affected more. 
Notice that the NS solution is much cleaner in the smooth parts close to the boundary, and the TR 
one is brighter than the original close to the center. This shows that the NS solutions reconstructs 
the low frequency modes better. 
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Boundary distance map The exact initial condition 




Figure 10: Example 2 with the trapping speed C3. Case 2: T = 4Tq with 10% random noise, (a): 
the boundary distance map. (b): the exact initial condition, (c): the time reversal solution, (d): 
the Neumann series solution. 
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Figure 11: Example 2 with the trapping speed C2- T = ATq. (a): the boundary distance map. (b): 
the exact initial condition, (c): the time reversal solution, (d): the Neumann series solution. 
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The velocity The velocity 




Figure 12: Sound speed models, (a): a discontinuous piecewise sound speed C4. (b): a non- 
piecewise constant discontinuous sound speed C5. 

9 Numerical results: discontinuous sound speed 

We present here numerical examples with two discontinuous sound speeds illustrated in Figure 12. 
The first one, that we denote by C4, is equal to 0.8 in the square [-1, l] 2 , then jumps by about a 
factor of 2 from the interior to the exterior, and then jumps to 1. The second one, C5, has similar 
jumps but inside the square [—1, l] 2 is variable, equal to the speed (32). 

If the speed jumps by a factor of 2 when going out of a square, all rays hitting the boundary at 
an angle less than 60 degrees are completely reflected. Those rays that hit the boundary at an angle 
greater than 30 degrees generate a reflected ray hitting the boundary again without a transmitted 
component, etc. Therefore, all rays hitting the boundary of the smaller square at angles between 
30 and 60 degrees are completely trapped for all times. 

9.1 Example 3: The Shepp-Logan phantom, Figures 13—14 
9.1.1 Piecewise constant discontinuous speed C4 

Figure 13: The sound speed is C4 given by Figure 12(a). We have To ~ 1.50 with T = 4Tq. 
The artifacts in the TR image are quite strong and can be explained by the way that singularities 
propagate in this case. In Figure 1, for example, the ray that exits on top carries a fraction of 
the energy only. When we reverse the time, at the first contact with the outer boundary of the 
"skull", this ray will create a reflected one (not seen in the figure) together with the transmitted 
one, shown there. The reflected one is not part of the actual graph that we are trying to invert. It 
will reflect off d£l, and go back to the interior of O, creating more artificial rays, etc. If we had an 
infinite time T, then those artificial rays will be canceled by other such artificial rays, leading to an 
exact reconstruction, as T — > 00 (at a very slow logarithmic rate for smooth /); this explains the 
artifacts in the TR image. The NS series expansion creates very few artifacts, and the few that are 
seen are mostly due to singularities near or on the "skull" . 
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Boundary distance map The exact initial condition 




Figure 13: Example 3 with the discontinuous sound speed C4. T = 4Tq. (a): the boundary distance 
map. (b): the exact initial condition, (c): the time reversal solution, (d): the Neumann series 
solution. 
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Boundary distance map The exact initial condition 




Figure 14: Example 3 with the discontinuous sound speed C5. (a): the boundary distance map. 
(b): the exact initial condition, (c): the time reversal solution, (d): the Neumann series solution. 

9.1.2 Discontinuous speed C5 

Figure 14: The sound speed C5 is given by Figure 12(b). Here, To sa 1.36, T = 4To. The speed 
is not so symmetric anymore and the artifacts in the TR image are still strong but more random. 
The variable speed inside improves the NS image. 

9.2 Example 4: Zebras, Figure 15 

Figure 15: The sound speed C5 is given by Figure 12(b). Here, To « 1.36, and we take T = 4Tq. 
As in Figure 14, the TR reconstruction (error 21.83%) contains many wrong singularities, while 
the NS image (error 8%, k = 16) is very clean. 

10 Numerical results: partial data, Figures 16—19 

We present here numerical examples with data given on three or two sides only. In the first case, 
we remove data on the right hand side of the square; in the second one, we remove data on the right 
hand side and on the bottom of the square. We also use a smooth cutoff. We present examples with 
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Figure 15: Example 4 with the discontinuous speed C5: Figure 12(b). (a): the boundary distance 
map. (b): the exact initial condition, (c): the time reversal solution, (d): the Neumann series 
solution. 



Qian, Stefanov, Uhlmann and Zhao 



Thermo- and Photo-acoustic Tomography 



30 



Boundary distance map The Neumann series solution 




Figure 16: Non-trapping speed c\ with partial data on two adjacent sides. T = 4.7. (a): the 
boundary distance map, data on two sides, (b): The Shepp-Logan phantom reconstruction. 

the non-trapping and the trapping speeds that we used before. Note that the notion of trapping 
changes with partial data since the stability depends whether all singularities can reach T for times 
< T. Still, removing the data on some of the sides in addition to the much longer times that signals 
need to reach T produces worse reconstructions when the speed is trapping. 

10.1 Partial data, non-trapping speed c l5 Figure 16 

The speed c\ is given by equation (32). In all three cases, T = 4.7. 

10.1.1 The Shepp-Logan phantom, Figure 16 (a), (b) 

We estimate To to be To ~ 2.5. It is greater than before because T is not the whole dO, anymore. 
We have data on two adjacent sides. If the speed were constant, the singularities below the lower- 
left-to-upper-right diagonal that would exit in the sides without measurements are invisible. That 
would affect jumps across surfaces in that triangle with normals parallel or nearly parallel to that 
diagonal, roughly speaking. The speed is variable but not too far from constant, and we see the 
expected behavior. Note that this is a unstable case, and our analysis does not exclude an even 
exponential divergence (in the low frequency part) but the error gets smaller up to k = 8 when the 
computation is stopped. 

10.1.2 Zebras, Figure 17 (a), (b) 

We have data on three sides. In this case, T\ ~ 1.35, and we would have stability if the speed were 
constant. However, the speed is not constant but in some sense, not too far from a constant one. 
There might be a small set of geodesies that enter through the right-hand side and exit there as 
well, thus creating instability. The reconstruction is very good, however, with an error of 6.11%, 
k = 10. The chosen T is slightly greater than what would be the stability time, but the result is 
computed without the contribution of those geodesies. 
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Boundary distance map The Neumann series solution 




Figure 17: A modified example with the non-trapping speed c\, partial data. T = 4.7. (a): the 
boundary distance map, data on three sides, (b): the reconstructed "zebras" image, (c): the 
boundary distance map, data on two sides, (d): the reconstructed "zebras" image. 

10.1.3 Zebras, Figure 17 (c), (d) 

We use data on two sides. This is a very unstable case and one can see artifacts where the invisible 
singularities lie — in the lower right-hand triangle, with slopes close to 1, roughly speaking (jumps 
across curves with slopes in a neighborhood of —1). The error is 10.6%, and k = 10. 

10.1.4 Zebras, Figure 18 (a), (b) 

The zebras image reconstructed in Figure 17 (d) does not have so many invisible singularities in 
this particular setup (data on two sides), however. For this reason, we present another example, 
Figure 17 (a)(b), with a modified image that shows the expected behavior of the visible and invisible 
singularities. 

10.2 Partial data, trapping sound speeds c 2 and c 2 , Figure 19 

The sound speed is C3 (the first two), and c 2 (the third example). 
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The exact initial condition The Neumann series solution 




Figure 18: Non-trapping speed c\ examples, partial data. T = 4.7. (a): the exact initial condition, 
(b): the reconstructed "zebras" image. 

10.2.1 Figure 19 (a), (b) 

Here, T = 4.93 > To ~ 1.3, with data on three sides. The "chaotic" trapping speed C3 makes the 
reconstruction worse than before. This is an unstable case because the trapping speed leaves many 
singularities invisible. As expected, the worst part is near the side with no observations due to 
geodesies that enter and exit through that side. There are invisible singularities everywhere, as 
well, due to the speed. 

10.2.2 Figure 19 (c), (d) 

Here, T = 4.93 > To « 2.1, with data on two sides, the same speed as above. As expected, the 
reconstruction is quite bad near the sides with no data. 

10.2.3 Figure 19 (e), (f) 

Here, T = 8.61 > To 2.6 with data on two sides, and the speed is c<i- The time T is larger than 
above but To is slightly larger as well. The reconstruction is better due to the larger time and 
(probably) due to the fact that the trapping region of this speed is farther away from the sides 
where no observations are done. Experiments with times T closer to that in the two examples 
above, not shown, still yield a better reconstruction with this speed. 

11 Conclusion 

We present new algorithms for reconstructing an unknown source in TAT and PAT based on the 
recent advances in understanding the theoretical nature of the problem. We work with variable 
sound speeds that might be also discontinuous across some surface. The latter problem arises in 
brain imaging. The new algorithm is based on an explicit formula in the form of a Neumann series. 
We present numerical examples with non-trapping, trapping and piecewise smooth speeds, as well 
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Boundary distance map 
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Figure 19: Examples with the trapping speeds C3 
boundary distance map, data on three sides, (b): 
the boundary distance map, data on two sides. 
k = 16. (e): the boundary distance map, data on 
T = 8.62, k= 16. 
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(the first two rows), and C2 (the last row), (a): the 
the reconstructed "zebras" image, T = 4.92. (c): 

(d): the reconstructed "zebras" image, T = 4.92, 
two sides, (f): the reconstructed "zebras" image, 
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as examples with data on a part of the boundary. These numerical examples demonstrate the 
robust performance of the new algorithm. 
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